##############################################################
#### 0. Load Packages & Read Data
##############################################################

rm(list=ls())

Sys.setlocale("LC_CTYPE", "en_US.UTF-8")

packages <- c("plyr", "stargazer", "dplyr", "tidyr", "nnet", "xtable", "dotwhisker", "ggplot2", "sjlabelled", "tidyverse", "labelled", "haven", "dplyr", "stargazer","cregg", "tidyverse","estimatr","clubSandwich","lmtest", "dotwhisker", "fastDummies", "cjoint", "monomvn", "stringr","survey", "estimatr", "coefplot", "glmnet", "haven")

package_installed <-
  sapply(packages, function(pack)
    pack %in% rownames(installed.packages()))

if (any(!package_installed)) {
  sapply(packages[!package_installed], install.packages)
}

sapply(packages, require, character.only = TRUE)

rm(packages,package_installed)

dat_conjoint_all <- readRDS("dat_conjoint_covid.rds")
names(dat_conjoint_all) <- gsub("_eng", "", names(dat_conjoint_all))

##############################################################
#### Figure 1: Effects of Immigrant Attributes on Probability of Supporting Admission
##############################################################

dat_conjoint_acme <- dat_conjoint_all[c("ResponseId", "selected", "residence", "trade", "regime", "covid.cases", "vaccine", "education", "korean", "occupation")]

acme <- cj(
  formula = selected ~ residence + trade + regime + covid.cases + vaccine + education + occupation +  korean,
  data = dat_conjoint_acme, 
  id = ~ ResponseId)

acme$feature <- factor(acme$feature, levels = c("regime", "trade", "residence", "covid.cases", "vaccine", "korean", "education", "occupation"))
acme$feature_label <- acme$feature
acme$feature_label <- gsub("residence", "Residence", acme$feature_label)
acme$feature_label <- gsub("trade", "Trade with South Korea", acme$feature_label)
acme$feature_label <- gsub("regime", "Regime type", acme$feature_label)
acme$feature_label <- gsub("covid.cases", "COVID-19 cases", acme$feature_label)
acme$feature_label <- gsub("vaccine", "COVID-19 vaccine", acme$feature_label)
acme$feature_label <- gsub("korean", "Language proficiency", acme$feature_label)
acme$feature_label <- gsub("education", "Education", acme$feature_label)
acme$feature_label <- gsub("occupation", "Employed sector", acme$feature_label)

acme$feature_label <- factor(acme$feature_label, levels = c("Employed sector", "Education", "Language proficiency", "COVID-19 vaccine", "COVID-19 cases", "Residence", "Trade with South Korea",  "Regime type"))

acme$level <- factor(acme$level, levels = c("Non-democracy", "Democracy", "Low volume of trade", "High volumne of trade", "Local community", "Expat community", "COVID-19: 5 per 10,000", "COVID-19: 50 per 10,000", "COVID-19: 500 per 10,000", "Unvaccinated", "Vaccinated", "No communication skills", "Basic communication skills", "Fluent communitation skills", "Middle school or below", "High school", "College or above", "Agriculture & Fishing", "Education", "Construction", "Finance", "Hospitality", "Manufacturing", "Science & Technology"))

pdf("Figure1.pdf")

acme %>% 
  ggplot(aes(x=level, y=estimate, color=feature_label)) +
  geom_point(position=position_dodge(width = .2)) +
  geom_linerange(aes(ymin=lower, ymax=upper), position=position_dodge(width = .2)) +
  geom_hline(yintercept=0, linetype="dashed", color="grey") +
  coord_flip() + xlab("Attribute levels") + ylab("Estimate") +
  #theme(legend.title = element_blank())  + 
  labs(color='Immigrant attributes')

dev.off()

##############################################################
#### Figure 2: Effects of Immigrant Attributes, By College Education
##############################################################

dat_conjoint_acme_educ <- dat_conjoint_all[c("ResponseId", "selected", "residence", "trade", "regime", "covid.cases", "vaccine", "education", "korean", "occupation", "college", "grad")]
names(dat_conjoint_acme_educ) <- gsub("", "", names(dat_conjoint_acme_educ))
dat_conjoint_acme_college <- subset(dat_conjoint_acme_educ, college == 1 | grad == 1)
dat_conjoint_acme_nocollege <- subset(dat_conjoint_acme_educ, college == 0 & grad == 0)

acme_college <- cj(
  formula = selected ~ residence + trade + regime + covid.cases + vaccine + education + occupation +  korean,
  data = dat_conjoint_acme_college, 
  id = ~ ResponseId)

acme_college$feature <- factor(acme_college$feature, levels = c("regime", "trade", "residence", "covid.cases", "vaccine", "korean", "education", "occupation"))
acme_college$feature_label <- acme_college$feature
acme_college$feature_label <- gsub("residence", "Residence", acme_college$feature_label)
acme_college$feature_label <- gsub("trade", "Trade with South Korea", acme_college$feature_label)
acme_college$feature_label <- gsub("regime", "Regime type", acme_college$feature_label)
acme_college$feature_label <- gsub("covid.cases", "COVID-19 cases", acme_college$feature_label)
acme_college$feature_label <- gsub("vaccine", "COVID-19 vaccine", acme_college$feature_label)
acme_college$feature_label <- gsub("korean", "Language proficiency", acme_college$feature_label)
acme_college$feature_label <- gsub("education", "Education", acme_college$feature_label)
acme_college$feature_label <- gsub("occupation", "Employed sector", acme_college$feature_label)

acme_college$feature_label <- factor(acme_college$feature_label, levels = c("Employed sector", "Education", "Language proficiency", "COVID-19 vaccine", "COVID-19 cases", "Residence", "Trade with South Korea",  "Regime type"))

acme_college$level <- factor(acme_college$level, levels = c("Non-democracy", "Democracy", "Low volume of trade", "High volumne of trade", "Local community", "Expat community", "COVID-19: 5 per 10,000", "COVID-19: 50 per 10,000", "COVID-19: 500 per 10,000", "Unvaccinated", "Vaccinated", "No communication skills", "Basic communication skills", "Fluent communitation skills", "Middle school or below", "High school", "College or above", "Agriculture & Fishing", "Education", "Construction", "Finance", "Hospitality", "Manufacturing", "Science & Technology"))
acme_college$College <- "College"

acme_nocollege <- cj(
  formula = selected ~ residence + trade + regime + covid.cases + vaccine + education + occupation +  korean,
  data = dat_conjoint_acme_nocollege, 
  id = ~ ResponseId)

acme_nocollege$feature <- factor(acme_nocollege$feature, levels = c("regime", "trade", "residence", "covid.cases", "vaccine", "korean", "education", "occupation"))
acme_nocollege$feature_label <- acme_nocollege$feature
acme_nocollege$feature_label <- gsub("residence", "Residence", acme_nocollege$feature_label)
acme_nocollege$feature_label <- gsub("trade", "Trade with South Korea", acme_nocollege$feature_label)
acme_nocollege$feature_label <- gsub("regime", "Regime type", acme_nocollege$feature_label)
acme_nocollege$feature_label <- gsub("covid.cases", "COVID-19 cases", acme_nocollege$feature_label)
acme_nocollege$feature_label <- gsub("vaccine", "COVID-19 vaccine", acme_nocollege$feature_label)
acme_nocollege$feature_label <- gsub("korean", "Language proficiency", acme_nocollege$feature_label)
acme_nocollege$feature_label <- gsub("education", "Education", acme_nocollege$feature_label)
acme_nocollege$feature_label <- gsub("occupation", "Employed sector", acme_nocollege$feature_label)

acme_nocollege$feature_label <- factor(acme_nocollege$feature_label, levels = c("Employed sector", "Education", "Language proficiency", "COVID-19 vaccine", "COVID-19 cases", "Residence", "Trade with South Korea",  "Regime type"))

acme_nocollege$level <- factor(acme_nocollege$level, levels = c("Non-democracy", "Democracy", "Low volume of trade", "High volumne of trade", "Local community", "Expat community", "COVID-19: 5 per 10,000", "COVID-19: 50 per 10,000", "COVID-19: 500 per 10,000", "Unvaccinated", "Vaccinated", "No communication skills", "Basic communication skills", "Fluent communitation skills", "Middle school or below", "High school", "College or above", "Agriculture & Fishing", "Education", "Construction", "Finance", "Hospitality", "Manufacturing", "Science & Technology"))
acme_nocollege$College <- "No College"

acme_educ <- rbind(acme_nocollege, acme_college)

pdf("Figure2.pdf")

acme_educ %>% 
  ggplot(aes(x=level, y=estimate, color=feature_label, shape=College)) +
  geom_point(position=position_dodge(width = .2)) +
  geom_linerange(aes(ymin=lower, ymax=upper), position=position_dodge(width = .2)) +
  geom_hline(yintercept=0, linetype="dashed", color="grey") +
  coord_flip() + xlab("Attribute levels") + ylab("Estimate") +
  #theme(legend.title = element_blank())  + 
  labs(color='Immigrant attributes', shape='Education')

dev.off()

##############################################################
#### Figure 3: Effects of Immigrant Attributes Conditional on Vaccination Status
##############################################################

dat_conjoint <- dat_conjoint_all[c("ResponseId", "selected", "residence", "trade", "regime", "covid.cases", "vaccine", "education", "korean", "occupation", "taskNum")]

dat_conjoint$task_pair <- paste0(as.character(dat_conjoint$ResponseId), as.character(dat_conjoint$taskNum), sep="_")

ANOVAfit <- CausalANOVA(formula= selected ~ residence + trade + regime + covid.cases + vaccine + education + occupation +  korean,
                        int2.formula = ~ residence:vaccine + trade:vaccine + regime:vaccine + covid.cases:vaccine + education:vaccine + occupation:vaccine + korean:vaccine,
                        data=dat_conjoint, 
                        pair.id=dat_conjoint$task_pair,
                        diff=TRUE,
                        cluster=dat_conjoint$ResponseId, nway=2)

gen_cond_eff<- function(fit,treat.fac,cond.fac){
  tfac<-treat.fac
  cfac<-cond.fac
  fit<-fit
  CE<-ConditionalEffect(fit, treat.fac=tfac, cond.fac=cfac)
  CE.plot <- CE$Conditional
  cond.level <- CE$cond.level
  treat.level <- CE$treat.level
  font <- CE.value <- CE.up <- CE.low <- Name<- c()
  #Name <- list()
  for(z in 1:length(CE.plot)){
    Name <- c(Name, c(paste(cond.fac, "=", cond.level[z],":",treat.level)))
    CE.value <- c(CE.value, c(CE.plot[[z]][,1]))
    CE.up <- c(CE.up, c(CE.plot[[z]][,4]))
    CE.low <- c(CE.low, c(CE.plot[[z]][,3]))
  }
  xlim.min <- min(na.omit(CE.low))
  xlim.max <- max(na.omit(CE.up))
  a<-data.frame(Name,CE.value,CE.up,CE.low) 
  return(a)
}

residence_reg<-gen_cond_eff(ANOVAfit,"residence", "vaccine")
trade_reg<-gen_cond_eff(ANOVAfit, "trade", "vaccine")
regime_reg<-gen_cond_eff(ANOVAfit,"regime", "vaccine")
covid_reg<-gen_cond_eff(ANOVAfit,"covid.cases", "vaccine")
education_reg<-gen_cond_eff(ANOVAfit,"education","vaccine")
occupation_reg<-gen_cond_eff(ANOVAfit,"occupation", "vaccine")
korean_reg<-gen_cond_eff(ANOVAfit,"korean", "vaccine")

causal_hetero<-rbind(residence_reg, trade_reg, regime_reg, covid_reg,education_reg, occupation_reg, korean_reg)
causal_hetero <-causal_hetero%>% mutate(Name1=Name) %>% tidyr::separate(Name1, c("type", "var"), " : ")

causal_hetero$var <- factor(causal_hetero$var, levels = c("Non-democracy", "Democracy", "Low volume of trade", "High volumne of trade", "Local community", "Expat community", "COVID-19: 5 per 10,000", "COVID-19: 50 per 10,000", "COVID-19: 500 per 10,000", "Unvaccinated", "Vaccinated", "No communication skills", "Basic communication skills", "Fluent communitation skills", "Middle school or below", "High school", "College or above", "Agriculture & Fishing", "Education", "Construction", "Finance", "Hospitality", "Manufacturing", "Science & Technology"))

causal_hetero$feature <- causal_hetero$var
causal_hetero$feature <- gsub('Non-democracy|Democracy', 'Regime type', causal_hetero$feature)
causal_hetero$feature <- gsub('Low volume of trade|High volumne of trade', 'Trade with South Korea', causal_hetero$feature)
causal_hetero$feature <- gsub('Local community|Expat community', 'Residence', causal_hetero$feature)
causal_hetero$feature <- gsub('COVID-19: 5 per 10,000|COVID-19: 50 per 10,000|COVID-19: 500 per 10,000', 'COVID-19 cases', causal_hetero$feature)
causal_hetero$feature <- gsub('No communication skills|Basic communication skills|Fluent communitation skills', 'Language proficiency', causal_hetero$feature)
causal_hetero$feature <- gsub('Agriculture & Fishing|Education|Construction|Finance|Hospitality|Manufacturing|Science & Technology', 'Employed sector', causal_hetero$feature)
causal_hetero$feature <- gsub('Middle school or below|High school|College or above', 'Education', causal_hetero$feature)

causal_hetero$feature <- factor(causal_hetero$feature, levels = c("Employed sector", "Education", "Language proficiency","COVID-19 cases", "Residence", "Trade with South Korea",  "Regime type"))
causal_hetero$type <- gsub("vaccine = ", "", causal_hetero$type)
causal_hetero$type <- factor(causal_hetero$type, levels = c("Vaccinated", "Unvaccinated"))

pdf("Figure3.pdf")

causal_hetero %>% 
  ggplot(aes(x=var, y=CE.value, color=as.factor(feature), shape=as.factor(type))) +
  geom_point(position=position_dodge(width = .2)) +
  geom_linerange(aes(ymin=CE.low, ymax=CE.up), position=position_dodge(width = .2)) +
  geom_hline(yintercept=0, linetype="dashed", color="grey") +
  coord_flip() + xlab("Attribute levels") + ylab("Estimate") +
  #theme(legend.title = element_blank())  + 
  labs(color='Immigrant attributes', shape="COVID-19 vaccine")

dev.off()

##############################################################
#### Figure A3: Robustness Test: An Alternative Outcome Variable (AMCE)
##############################################################

dat_conjoint_quarantine <- dat_conjoint_all[c("ResponseId", "selected", "residence", "trade", "regime", "covid.cases", "vaccine", "education", "korean", "occupation", "restriction")]
names(dat_conjoint_quarantine) <- gsub("", "", names(dat_conjoint_quarantine))

dat_conjoint_quarantine$entry_allowed <- ifelse(dat_conjoint_quarantine$restriction=="입국 허용", 1, 0)
dat_conjoint_quarantine$quarantine <- ifelse(dat_conjoint_quarantine$restriction=="입국 금지", 0, 1)

acme_quarantine <- cj(
  formula = quarantine ~ residence + trade + regime + covid.cases + vaccine + education + occupation +  korean,
  data = dat_conjoint_quarantine, 
  id = ~ ResponseId)

acme_quarantine$feature <- factor(acme_quarantine$feature, levels = c("regime", "trade", "residence", "covid.cases", "vaccine", "korean", "education", "occupation"))
acme_quarantine$feature_label <- acme_quarantine$feature
acme_quarantine$feature_label <- gsub("residence", "Residence", acme_quarantine$feature_label)
acme_quarantine$feature_label <- gsub("trade", "Trade with South Korea", acme_quarantine$feature_label)
acme_quarantine$feature_label <- gsub("regime", "Regime type", acme_quarantine$feature_label)
acme_quarantine$feature_label <- gsub("covid.cases", "COVID-19 cases", acme_quarantine$feature_label)
acme_quarantine$feature_label <- gsub("vaccine", "COVID-19 vaccine", acme_quarantine$feature_label)
acme_quarantine$feature_label <- gsub("korean", "Language proficiency", acme_quarantine$feature_label)
acme_quarantine$feature_label <- gsub("education", "Education", acme_quarantine$feature_label)
acme_quarantine$feature_label <- gsub("occupation", "Employed sector", acme_quarantine$feature_label)

acme_quarantine$feature_label <- factor(acme_quarantine$feature_label, levels = c("Employed sector", "Education", "Language proficiency", "COVID-19 vaccine", "COVID-19 cases", "Residence", "Trade with South Korea",  "Regime type"))

acme_quarantine$level <- factor(acme_quarantine$level, levels = c("Non-democracy", "Democracy", "Low volume of trade", "High volumne of trade", "Local community", "Expat community", "COVID-19: 5 per 10,000", "COVID-19: 50 per 10,000", "COVID-19: 500 per 10,000", "Unvaccinated", "Vaccinated", "No communication skills", "Basic communication skills", "Fluent communitation skills", "Middle school or below", "High school", "College or above", "Agriculture & Fishing", "Education", "Construction", "Finance", "Hospitality", "Manufacturing", "Science & Technology"))

pdf("FigureA3.pdf")

acme_quarantine %>% 
  ggplot(aes(x=level, y=estimate, color=feature_label)) +
  geom_point(position=position_dodge(width = .2)) +
  geom_linerange(aes(ymin=lower, ymax=upper), position=position_dodge(width = .2)) +
  geom_hline(yintercept=0, linetype="dashed", color="grey") +
  coord_flip() + xlab("Attribute levels") + ylab("Estimate") +
  #theme(legend.title = element_blank())  + 
  labs(color='Immigrant attributes')

dev.off()


##############################################################
#### Figure A4: Robustness Test: An Alternative Outcome Variable (Multinomial Logit)
##############################################################

dat_conjoint_multi <- dat_conjoint_all
dat_conjoint_multi <- fastDummies::dummy_cols(dat_conjoint_multi, select_columns = c("residence", "occupation", "trade", "regime", "covid.cases", "vaccine", "education", "korean"))

names(dat_conjoint_multi) <- gsub("[a-z\\.]+_", "", names(dat_conjoint_multi))
dat_conjoint_multi$restriction_type <- relevel(factor(dat_conjoint_multi$restriction), ref = 3)

data = dat_conjoint_multi
formula =  restriction_type ~ `Expat community` + `Construction` + `Education` + `Finance` + `Hospitality` + `Manufacturing` + `Science & Technology` + `High volumne of trade` + `Democracy` +`COVID-19: 50 per 10,000` + `COVID-19: 500 per 10,000` + `Vaccinated` + `College or above` + `High school` + `Basic communication skills` + `Fluent communitation skills` 

test <- multinom(formula, data = data)
df_mnl <- tidy(test) %>% mutate(model = y.level)  %>% filter(term !='(Intercept)')
df_mnl$upper <- df_mnl$estimate+df_mnl$std.error*2
df_mnl$lower <- df_mnl$estimate-df_mnl$std.error*2

df_mnl$feature <- c("Residence", "Employed sector", "Employed sector", "Employed sector", "Employed sector", "Employed sector", "Employed sector", "Trade with South Korea", "Regime type", "COVID-19 cases", "COVID-19 cases", "COVID-19 vaccine", "Education", "Education", "Language proficiency", "Language proficiency",
                    "Residence", "Employed sector", "Employed sector", "Employed sector", "Employed sector", "Employed sector", "Employed sector", "Trade with South Korea", "Regime type", "COVID-19 cases", "COVID-19 cases", "COVID-19 vaccine", "Education", "Education", "Language proficiency", "Language proficiency" )

df_mnl$feature <- factor(df_mnl$feature, levels = c("Employed sector", "Education", "Language proficiency", "COVID-19 vaccine", "COVID-19 cases", "Residence", "Trade with South Korea",  "Regime type"))
df_mnl$term <- gsub("`", "", df_mnl$term )

df_mnl$Choice <- c(rep("Entry without quarantine", 16), rep("Entry after 14 days of quarantine", 16))

df_mnl$term <- factor(df_mnl$term, levels = c("Democracy", "High volumne of trade", "Expat community", "COVID-19: 50 per 10,000", "COVID-19: 500 per 10,000", "Vaccinated", "Basic communication skills", "Fluent communitation skills", "High school", "College or above", "Education", "Construction", "Finance", "Hospitality", "Manufacturing", "Science & Technology"))

pdf("FigureA4.pdf")

df_mnl %>% 
  ggplot(aes(x=term, y=estimate, color=feature, shape = Choice)) +
  geom_point(position=position_dodge(width = .2)) +
  geom_linerange(aes(ymin=lower, ymax=upper), position=position_dodge(width = .2)) +
  geom_hline(yintercept=0, linetype="dashed", color="grey") +
  coord_flip() + xlab("Attribute levels") + ylab("Estimate") +
  #theme(legend.title = element_blank())  + 
  labs(color='Immigrant attributes')

dev.off()

##############################################################
#### Figure A5: Robustness Test: Excluding Inattentive Responses (1)
##############################################################

dat_conjoint_acme_robust <- subset(dat_conjoint_all, Duration__in_seconds_ >= 300)

dat_conjoint_acme_robust$Duration__in_seconds_ <- NULL
dat_conjoint_acme_robust$same_answer <- NULL

acme_robust1 <- cj(
  formula = selected ~ residence + trade + regime + covid.cases + vaccine + education + occupation +  korean,
  data = dat_conjoint_acme_robust, 
  id = ~ ResponseId)

acme_robust1$feature <- factor(acme_robust1$feature, levels = c("regime", "trade", "residence", "covid.cases", "vaccine", "korean", "education", "occupation"))
acme_robust1$feature_label <- acme_robust1$feature
acme_robust1$feature_label <- gsub("residence", "Residence", acme_robust1$feature_label)
acme_robust1$feature_label <- gsub("trade", "Trade with South Korea", acme_robust1$feature_label)
acme_robust1$feature_label <- gsub("regime", "Regime type", acme_robust1$feature_label)
acme_robust1$feature_label <- gsub("covid.cases", "COVID-19 cases", acme_robust1$feature_label)
acme_robust1$feature_label <- gsub("vaccine", "COVID-19 vaccine", acme_robust1$feature_label)
acme_robust1$feature_label <- gsub("korean", "Language proficiency", acme_robust1$feature_label)
acme_robust1$feature_label <- gsub("education", "Education", acme_robust1$feature_label)
acme_robust1$feature_label <- gsub("occupation", "Employed sector", acme_robust1$feature_label)

acme_robust1$feature_label <- factor(acme_robust1$feature_label, levels = c("Employed sector", "Education", "Language proficiency", "COVID-19 vaccine", "COVID-19 cases", "Residence", "Trade with South Korea",  "Regime type"))

acme_robust1$level <- factor(acme_robust1$level, levels = c("Non-democracy", "Democracy", "Low volume of trade", "High volumne of trade", "Local community", "Expat community", "COVID-19: 5 per 10,000", "COVID-19: 50 per 10,000", "COVID-19: 500 per 10,000", "Unvaccinated", "Vaccinated", "No communication skills", "Basic communication skills", "Fluent communitation skills", "Middle school or below", "High school", "College or above", "Agriculture & Fishing", "Education", "Construction", "Finance", "Hospitality", "Manufacturing", "Science & Technology"))

pdf("FigureA5.pdf")

acme_robust1 %>% 
  ggplot(aes(x=level, y=estimate, color=feature_label)) +
  geom_point(position=position_dodge(width = .2)) +
  geom_linerange(aes(ymin=lower, ymax=upper), position=position_dodge(width = .2)) +
  geom_hline(yintercept=0, linetype="dashed", color="grey") +
  coord_flip() + xlab("Attribute levels") + ylab("Estimate") +
  #theme(legend.title = element_blank())  + 
  labs(color='Immigrant attributes')

dev.off()

##############################################################
#### Figure A6: Robustness Test: Excluding Inattentive Responses (2)
##############################################################

dat_conjoint_acme_robust <- subset(dat_conjoint_all, same_answer != 1)
dat_conjoint_acme_robust$Duration__in_seconds_ <- NULL
dat_conjoint_acme_robust$same_answer <- NULL

acme_robust2 <- cj(
  formula = selected ~ residence + trade + regime + covid.cases + vaccine + education + occupation +  korean,
  data = dat_conjoint_acme_robust, 
  id = ~ ResponseId)

acme_robust2$feature <- factor(acme_robust2$feature, levels = c("regime", "trade", "residence", "covid.cases", "vaccine", "korean", "education", "occupation"))
acme_robust2$feature_label <- acme_robust2$feature
acme_robust2$feature_label <- gsub("residence", "Residence", acme_robust2$feature_label)
acme_robust2$feature_label <- gsub("trade", "Trade with South Korea", acme_robust2$feature_label)
acme_robust2$feature_label <- gsub("regime", "Regime type", acme_robust2$feature_label)
acme_robust2$feature_label <- gsub("covid.cases", "COVID-19 cases", acme_robust2$feature_label)
acme_robust2$feature_label <- gsub("vaccine", "COVID-19 vaccine", acme_robust2$feature_label)
acme_robust2$feature_label <- gsub("korean", "Language proficiency", acme_robust2$feature_label)
acme_robust2$feature_label <- gsub("education", "Education", acme_robust2$feature_label)
acme_robust2$feature_label <- gsub("occupation", "Employed sector", acme_robust2$feature_label)

acme_robust2$feature_label <- factor(acme_robust2$feature_label, levels = c("Employed sector", "Education", "Language proficiency", "COVID-19 vaccine", "COVID-19 cases", "Residence", "Trade with South Korea",  "Regime type"))

acme_robust2$level <- factor(acme_robust2$level, levels = c("Non-democracy", "Democracy", "Low volume of trade", "High volumne of trade", "Local community", "Expat community", "COVID-19: 5 per 10,000", "COVID-19: 50 per 10,000", "COVID-19: 500 per 10,000", "Unvaccinated", "Vaccinated", "No communication skills", "Basic communication skills", "Fluent communitation skills", "Middle school or below", "High school", "College or above", "Agriculture & Fishing", "Education", "Construction", "Finance", "Hospitality", "Manufacturing", "Science & Technology"))

pdf("FigureA6.pdf")

acme_robust2 %>% 
  ggplot(aes(x=level, y=estimate, color=feature_label)) +
  geom_point(position=position_dodge(width = .2)) +
  geom_linerange(aes(ymin=lower, ymax=upper), position=position_dodge(width = .2)) +
  geom_hline(yintercept=0, linetype="dashed", color="grey") +
  coord_flip() + xlab("Attribute levels") + ylab("Estimate") +
  #theme(legend.title = element_blank())  + 
  labs(color='Immigrant attributes')

dev.off()

##############################################################
#### Figure A7: Effects of Immigrant Attributes on Probability of Supporting Admission: Comparison of Estimates from Different Models
##############################################################

acme$model <- "Main"
acme_robust1$model <- "Robustness 1"
acme_robust2$model <- "Robustness 2"

acme_all <- rbind(acme, acme_robust1)
acme_all <- rbind(acme_all, acme_robust2)

pdf("FigureA7.pdf")

acme_all %>% 
  ggplot(aes(x=level, y=estimate, color=feature_label, shape=model)) +
  geom_point(position=position_dodge(width = .2), alpha = 0.7) +
  geom_linerange(aes(ymin=lower, ymax=upper), position=position_dodge(width = .2), alpha = 0.7) +
  geom_hline(yintercept=0, linetype="dashed", color="grey") +
  coord_flip() + xlab("Attribute levels") + ylab("Estimate") +
  #theme(legend.title = element_blank())  + 
  labs(color='Immigrant attributes', shape = 'Model')

dev.off()

##############################################################
#### Figure A8: Effects of Immigrant Attributes on Probability of Supporting Admission, By Gender
##############################################################

dat_conjoint_acme_female <- subset(dat_conjoint_all, female == 1)
dat_conjoint_acme_male <- subset(dat_conjoint_all, female == 0)

acme_female <- cj(
  formula = selected ~ residence + trade + regime + covid.cases + vaccine + education + occupation +  korean,
  data = dat_conjoint_acme_female, 
  id = ~ ResponseId)

acme_female$feature <- factor(acme_female$feature, levels = c("regime", "trade", "residence", "covid.cases", "vaccine", "korean", "education", "occupation"))
acme_female$feature_label <- acme_female$feature
acme_female$feature_label <- gsub("residence", "Residence", acme_female$feature_label)
acme_female$feature_label <- gsub("trade", "Trade with South Korea", acme_female$feature_label)
acme_female$feature_label <- gsub("regime", "Regime type", acme_female$feature_label)
acme_female$feature_label <- gsub("covid.cases", "COVID-19 cases", acme_female$feature_label)
acme_female$feature_label <- gsub("vaccine", "COVID-19 vaccine", acme_female$feature_label)
acme_female$feature_label <- gsub("korean", "Language proficiency", acme_female$feature_label)
acme_female$feature_label <- gsub("education", "Education", acme_female$feature_label)
acme_female$feature_label <- gsub("occupation", "Employed sector", acme_female$feature_label)

acme_female$feature_label <- factor(acme_female$feature_label, levels = c("Employed sector", "Education", "Language proficiency", "COVID-19 vaccine", "COVID-19 cases", "Residence", "Trade with South Korea",  "Regime type"))

acme_female$level <- factor(acme_female$level, levels = c("Non-democracy", "Democracy", "Low volume of trade", "High volumne of trade", "Local community", "Expat community", "COVID-19: 5 per 10,000", "COVID-19: 50 per 10,000", "COVID-19: 500 per 10,000", "Unvaccinated", "Vaccinated", "No communication skills", "Basic communication skills", "Fluent communitation skills", "Middle school or below", "High school", "College or above", "Agriculture & Fishing", "Education", "Construction", "Finance", "Hospitality", "Manufacturing", "Science & Technology"))
acme_female$Gender <- "Female"

acme_male <- cj(
  formula = selected ~ residence + trade + regime + covid.cases + vaccine + education + occupation +  korean,
  data = dat_conjoint_acme_male, 
  id = ~ ResponseId)

acme_male$feature <- factor(acme_male$feature, levels = c("regime", "trade", "residence", "covid.cases", "vaccine", "korean", "education", "occupation"))
acme_male$feature_label <- acme_male$feature
acme_male$feature_label <- gsub("residence", "Residence", acme_male$feature_label)
acme_male$feature_label <- gsub("trade", "Trade with South Korea", acme_male$feature_label)
acme_male$feature_label <- gsub("regime", "Regime type", acme_male$feature_label)
acme_male$feature_label <- gsub("covid.cases", "COVID-19 cases", acme_male$feature_label)
acme_male$feature_label <- gsub("vaccine", "COVID-19 vaccine", acme_male$feature_label)
acme_male$feature_label <- gsub("korean", "Language proficiency", acme_male$feature_label)
acme_male$feature_label <- gsub("education", "Education", acme_male$feature_label)
acme_male$feature_label <- gsub("occupation", "Employed sector", acme_male$feature_label)

acme_male$feature_label <- factor(acme_male$feature_label, levels = c("Employed sector", "Education", "Language proficiency", "COVID-19 vaccine", "COVID-19 cases", "Residence", "Trade with South Korea",  "Regime type"))

acme_male$level <- factor(acme_male$level, levels = c("Non-democracy", "Democracy", "Low volume of trade", "High volumne of trade", "Local community", "Expat community", "COVID-19: 5 per 10,000", "COVID-19: 50 per 10,000", "COVID-19: 500 per 10,000", "Unvaccinated", "Vaccinated", "No communication skills", "Basic communication skills", "Fluent communitation skills", "Middle school or below", "High school", "College or above", "Agriculture & Fishing", "Education", "Construction", "Finance", "Hospitality", "Manufacturing", "Science & Technology"))
acme_male$Gender <- "Male"

acme_gender <- rbind(acme_male, acme_female)

pdf("FigureA8.pdf")

acme_gender %>% 
  ggplot(aes(x=level, y=estimate, color=feature_label, shape=Gender)) +
  geom_point(position=position_dodge(width = .2)) +
  geom_linerange(aes(ymin=lower, ymax=upper), position=position_dodge(width = .2)) +
  geom_hline(yintercept=0, linetype="dashed", color="grey") +
  coord_flip() + xlab("Attribute levels") + ylab("Estimate") +
  #theme(legend.title = element_blank())  + 
  labs(color='Immigrant attributes', shape='Gender')

dev.off()

##############################################################
#### Figure A9: Effects of Immigrant Attributes on Probability of Supporting Admission, By Age Groups
##############################################################

dat_conjoint_acme_2030 <- subset(dat_conjoint_all, Age_Group == "18-29" | Age_Group == "30-39")
dat_conjoint_acme_4050 <- subset(dat_conjoint_all, Age_Group == "40-49" | Age_Group == "50-59")
dat_conjoint_acme_60 <- subset(dat_conjoint_all, Age_Group == "60-69")

acme_2030 <- cj(
  formula = selected ~ residence + trade + regime + covid.cases + vaccine + education + occupation +  korean,
  data = dat_conjoint_acme_2030, 
  id = ~ ResponseId)

acme_2030$feature <- factor(acme_2030$feature, levels = c("regime", "trade", "residence", "covid.cases", "vaccine", "korean", "education", "occupation"))
acme_2030$feature_label <- acme_2030$feature
acme_2030$feature_label <- gsub("residence", "Residence", acme_2030$feature_label)
acme_2030$feature_label <- gsub("trade", "Trade with South Korea", acme_2030$feature_label)
acme_2030$feature_label <- gsub("regime", "Regime type", acme_2030$feature_label)
acme_2030$feature_label <- gsub("covid.cases", "COVID-19 cases", acme_2030$feature_label)
acme_2030$feature_label <- gsub("vaccine", "COVID-19 vaccine", acme_2030$feature_label)
acme_2030$feature_label <- gsub("korean", "Language proficiency", acme_2030$feature_label)
acme_2030$feature_label <- gsub("education", "Education", acme_2030$feature_label)
acme_2030$feature_label <- gsub("occupation", "Employed sector", acme_2030$feature_label)

acme_2030$feature_label <- factor(acme_2030$feature_label, levels = c("Employed sector", "Education", "Language proficiency", "COVID-19 vaccine", "COVID-19 cases", "Residence", "Trade with South Korea",  "Regime type"))

acme_2030$level <- factor(acme_2030$level, levels = c("Non-democracy", "Democracy", "Low volume of trade", "High volumne of trade", "Local community", "Expat community", "COVID-19: 5 per 10,000", "COVID-19: 50 per 10,000", "COVID-19: 500 per 10,000", "Unvaccinated", "Vaccinated", "No communication skills", "Basic communication skills", "Fluent communitation skills", "Middle school or below", "High school", "College or above", "Agriculture & Fishing", "Education", "Construction", "Finance", "Hospitality", "Manufacturing", "Science & Technology"))
acme_2030$Age_Group <- "Aged 20-39"

acme_4050 <- cj(
  formula = selected ~ residence + trade + regime + covid.cases + vaccine + education + occupation +  korean,
  data = dat_conjoint_acme_4050, 
  id = ~ ResponseId)

acme_4050$feature <- factor(acme_4050$feature, levels = c("regime", "trade", "residence", "covid.cases", "vaccine", "korean", "education", "occupation"))
acme_4050$feature_label <- acme_4050$feature
acme_4050$feature_label <- gsub("residence", "Residence", acme_4050$feature_label)
acme_4050$feature_label <- gsub("trade", "Trade with South Korea", acme_4050$feature_label)
acme_4050$feature_label <- gsub("regime", "Regime type", acme_4050$feature_label)
acme_4050$feature_label <- gsub("covid.cases", "COVID-19 cases", acme_4050$feature_label)
acme_4050$feature_label <- gsub("vaccine", "COVID-19 vaccine", acme_4050$feature_label)
acme_4050$feature_label <- gsub("korean", "Language proficiency", acme_4050$feature_label)
acme_4050$feature_label <- gsub("education", "Education", acme_4050$feature_label)
acme_4050$feature_label <- gsub("occupation", "Employed sector", acme_4050$feature_label)

acme_4050$feature_label <- factor(acme_4050$feature_label, levels = c("Employed sector", "Education", "Language proficiency", "COVID-19 vaccine", "COVID-19 cases", "Residence", "Trade with South Korea",  "Regime type"))

acme_4050$level <- factor(acme_4050$level, levels = c("Non-democracy", "Democracy", "Low volume of trade", "High volumne of trade", "Local community", "Expat community", "COVID-19: 5 per 10,000", "COVID-19: 50 per 10,000", "COVID-19: 500 per 10,000", "Unvaccinated", "Vaccinated", "No communication skills", "Basic communication skills", "Fluent communitation skills", "Middle school or below", "High school", "College or above", "Agriculture & Fishing", "Education", "Construction", "Finance", "Hospitality", "Manufacturing", "Science & Technology"))
acme_4050$Age_Group <- "Aged 40-59"

acme_60 <- cj(
  formula = selected ~ residence + trade + regime + covid.cases + vaccine + education + occupation +  korean,
  data = dat_conjoint_acme_60, 
  id = ~ ResponseId)

acme_60$feature <- factor(acme_60$feature, levels = c("regime", "trade", "residence", "covid.cases", "vaccine", "korean", "education", "occupation"))
acme_60$feature_label <- acme_60$feature
acme_60$feature_label <- gsub("residence", "Residence", acme_60$feature_label)
acme_60$feature_label <- gsub("trade", "Trade with South Korea", acme_60$feature_label)
acme_60$feature_label <- gsub("regime", "Regime type", acme_60$feature_label)
acme_60$feature_label <- gsub("covid.cases", "COVID-19 cases", acme_60$feature_label)
acme_60$feature_label <- gsub("vaccine", "COVID-19 vaccine", acme_60$feature_label)
acme_60$feature_label <- gsub("korean", "Language proficiency", acme_60$feature_label)
acme_60$feature_label <- gsub("education", "Education", acme_60$feature_label)
acme_60$feature_label <- gsub("occupation", "Employed sector", acme_60$feature_label)

acme_60$feature_label <- factor(acme_60$feature_label, levels = c("Employed sector", "Education", "Language proficiency", "COVID-19 vaccine", "COVID-19 cases", "Residence", "Trade with South Korea",  "Regime type"))

acme_60$level <- factor(acme_60$level, levels = c("Non-democracy", "Democracy", "Low volume of trade", "High volumne of trade", "Local community", "Expat community", "COVID-19: 5 per 10,000", "COVID-19: 50 per 10,000", "COVID-19: 500 per 10,000", "Unvaccinated", "Vaccinated", "No communication skills", "Basic communication skills", "Fluent communitation skills", "Middle school or below", "High school", "College or above", "Agriculture & Fishing", "Education", "Construction", "Finance", "Hospitality", "Manufacturing", "Science & Technology"))
acme_60$Age_Group <- "Aged 60+"

acme_age <- rbind(acme_2030, acme_4050)
acme_age <- rbind(acme_age, acme_60)

pdf("FigureA9.pdf")

acme_age %>% 
  ggplot(aes(x=level, y=estimate, color=feature_label, shape=Age_Group)) +
  geom_point(position=position_dodge(width = .2)) +
  geom_linerange(aes(ymin=lower, ymax=upper), position=position_dodge(width = .2)) +
  geom_hline(yintercept=0, linetype="dashed", color="grey") +
  coord_flip() + xlab("Attribute levels") + ylab("Estimate") +
  #theme(legend.title = element_blank())  + 
  labs(color='Immigrant attributes', shape='Age Group')

dev.off()

